function F = sys_irs(x,e2,ell,theta,sigma,pi,rho,eta,i,M,A,B)
% partial equilibrium, given e2
q_star = 1;

%-------------------------------------------------------------------------
% core variables
%-------------------------------------------------------------------------
% segmented markets
q0A  = x(1);
q1A  = x(2);
q0B  = x(3);
q1Bn = x(4);
q1Bd = x(5);
eC   = x(6);
eNn  = x(7);
eNd  = x(8);
% consolidated market
q0   = x(9);
q1n  = x(10); 
q1d  = x(11); 
% d's
dA   = x(12);
dB   = x(13);
dA2  = x(14);
dB2  = x(15);
% phi
phi  = 1/M * ((1-e2)*(eC*q0A+(1-eC)*q0B)+e2*q0);
% liquidity premia
L_A  = x(16);
L_A2 = x(17);
L_B_rho_bar  = x(18);
L_B2_rho_bar = x(19);

%-------------------------------------------------------------------------
% matching probabilities
%-------------------------------------------------------------------------
% segmented markets
a_CA_n = (1-e2)^eta * eNn*(1-ell)/(eC*ell+eNn*(1-ell))^(1-eta);
a_CB_n = (1-e2)^eta * (1-eNn)*(1-ell)/((1-eC)*ell+(1-eNn)*(1-ell))^(1-eta);
a_NA_n = (1-e2)^eta * eC*ell/(eC*ell+eNn*(1-ell))^(1-eta);
a_NB_n = (1-e2)^eta * (1-eC)*ell/((1-eC)*ell+(1-eNn)*(1-ell))^(1-eta);
a_CA_d = (1-e2)^eta * eNd*(1-ell)/(eC*ell+eNd*(1-ell))^(1-eta);
a_CB_d = (1-e2)^eta * (1-eNd)*(1-ell)/((1-eC)*ell+(1-eNd)*(1-ell))^(1-eta);
a_NA_d = (1-e2)^eta * eC*ell/(eC*ell+eNd*(1-ell))^(1-eta);
a_NB_d = (1-e2)^eta * (1-eC)*ell/((1-eC)*ell+(1-eNd)*(1-ell))^(1-eta);
% consolidated market
a_C2   = e2^eta * (1-ell);
a_N2   = e2^eta * ell;

%-------------------------------------------------------------------------
% segmented markets
%-------------------------------------------------------------------------
F(1) = - i + ell*(1-(pi*a_CA_n+(1-pi)*a_CA_d)*theta/w(q1A,sigma,theta))*(mu(q0A,sigma)-1) ...
           + ell*   (pi*a_CA_n+(1-pi)*a_CA_d)*theta/w(q1A,sigma,theta) *(mu(q1A,sigma)-1);
F(2) = - i + ell*(1-pi*a_CB_n*theta/w(q1Bn,sigma,theta)-(1-pi)*a_CB_d*theta/w(q1Bd,sigma,theta))*(mu(q0B,sigma)-1) ...
           + ell*pi*a_CB_n*theta/w(q1Bn,sigma,theta)*(mu(q1Bn,sigma)-1) + ell*(1-pi)*a_CB_d*theta/w(q1Bd,sigma,theta)*(mu(q1Bd,sigma)-1);

F(3) = - q1A  + min(q_star, q0A + (phi*dA     - (1-theta)*(u(q1A, sigma)-u(q0A,sigma)))/theta);
F(4) = - q1Bn + min(q_star, q0B + (phi*dB     - (1-theta)*(u(q1Bn,sigma)-u(q0B,sigma)))/theta);
F(5) = - q1Bd + min(q_star, q0B + (phi*dB*rho - (1-theta)*(u(q1Bd,sigma)-u(q0B,sigma)))/theta);

F(6) = dA * (- L_A + ell*(pi*a_CA_n+(1-pi)*a_CA_d)*theta/w(q1A,sigma,theta)*(mu(q1A,sigma)-1));
F(7) = dB * (- L_B_rho_bar + ell*pi*a_CB_n*theta/w(q1Bn,sigma,theta)*(mu(q1Bn,sigma)-1) + ell*(1-pi)*rho*a_CB_d*theta/w(q1Bd,sigma,theta)*(mu(q1Bd,sigma)-1));
% L_A         = ell*(pi*a_CA_n+(1-pi)*a_CA_d)*theta/w(q1A,sigma,theta)*(mu(q1A,sigma)-1);
% L_B_rho_bar = ell*pi*a_CB_n*theta/w(q1Bn,sigma,theta)*(mu(q1Bn,sigma)-1) + ell*(1-pi)*rho*a_CB_d*theta/w(q1Bd,sigma,theta)*(mu(q1Bd,sigma)-1);
S_CA   = theta*(u(q1A, sigma)-u(q0A,sigma)-q1A +q0A);
S_CB_n = theta*(u(q1Bn,sigma)-u(q0B,sigma)-q1Bn+q0B);
S_CB_d = theta*(u(q1Bd,sigma)-u(q0B,sigma)-q1Bd+q0B);
S_CA_tilde = - i*q0A - L_A*phi*dA + ell*(u(q0A,sigma)-q0A+(pi*a_CA_n+(1-pi)*a_CA_d)*S_CA);
S_CB_tilde = - i*q0B - L_B_rho_bar*phi*dB + ell*(u(q0B,sigma)-q0B+pi*a_CB_n*S_CB_n+(1-pi)*a_CB_d*S_CB_d);

F(8) = S_CA_tilde - S_CB_tilde;

S_NA   = (1-theta)*(u(q1A, sigma)-u(q0A,sigma)-q1A +q0A);
S_NB_n = (1-theta)*(u(q1Bn,sigma)-u(q0B,sigma)-q1Bn+q0B);
S_NB_d = (1-theta)*(u(q1Bd,sigma)-u(q0B,sigma)-q1Bd+q0B);

F(9)  = a_NA_n*S_NA - a_NB_n*S_NB_n;
F(10) = a_NA_d*S_NA - a_NB_d*S_NB_d;

S1 = S_CA_tilde + (1-ell)*pi*a_NA_n*S_NA + (1-ell)*(1-pi)*a_NA_d*S_NA;
S1 = S_CB_tilde + (1-ell)*pi*a_NB_n*S_NB_n + (1-ell)*(1-pi)*a_NB_d*S_NB_d;

%-------------------------------------------------------------------------
% consolidated market
%-------------------------------------------------------------------------
F(11) = - i + ell*(1-pi*a_C2*theta/w(q1n,sigma,theta)-(1-pi)*a_C2*theta/w(q1d,sigma,theta))*(mu(q0,sigma)-1) ...
            + ell*pi*a_C2*theta/w(q1n,sigma,theta)*(mu(q1n,sigma)-1) + ell*(1-pi)*a_C2*theta/w(q1d,sigma,theta)*(mu(q1d,sigma)-1);

F(12) = - q1n + min(q_star, q0 + (phi*dA2 + phi*dB2     - (1-theta)*(u(q1n,sigma)-u(q0,sigma)))/theta);
F(13) = - q1d + min(q_star, q0 + (phi*dA2 + phi*dB2*rho - (1-theta)*(u(q1d,sigma)-u(q0,sigma)))/theta);

F(14) = dA2 * (- L_A2 + ell*pi*a_C2*theta/w(q1n,sigma,theta)*(mu(q1n,sigma)-1) + ell*(1-pi)*a_C2*theta/w(q1d,sigma,theta)*(mu(q1d,sigma)-1));
F(15) = dB2 * (- L_B2_rho_bar + ell*pi*a_C2*theta/w(q1n,sigma,theta)*(mu(q1n,sigma)-1) + ell*(1-pi)*rho*a_C2*theta/w(q1d,sigma,theta)*(mu(q1d,sigma)-1));
% L_A2         = ell*pi*a_C2*theta/w(q1n,sigma,theta)*(mu(q1n,sigma)-1) + ell*(1-pi)*a_C2*theta/w(q1d,sigma,theta)*(mu(q1d,sigma)-1);
% L_B2_rho_bar = ell*pi*a_C2*theta/w(q1n,sigma,theta)*(mu(q1n,sigma)-1) + ell*(1-pi)*rho*a_C2*theta/w(q1d,sigma,theta)*(mu(q1d,sigma)-1);
S_C2_n = theta*(u(q1n,sigma)-u(q0,sigma)-q1n+q0);
S_C2_d = theta*(u(q1d,sigma)-u(q0,sigma)-q1d+q0);
S_C2_tilde = - i*q0 - L_A2*phi*dA2 - L_B2_rho_bar*phi*dB2 + ell*(u(q0,sigma)-q0+pi*a_C2*S_C2_n+(1-pi)*a_C2*S_C2_d);

S_N2_n = (1-theta)*(u(q1n,sigma)-u(q0,sigma)-q1n+q0);
S_N2_d = (1-theta)*(u(q1d,sigma)-u(q0,sigma)-q1d+q0);

S2 = S_C2_tilde + (1-ell)*a_N2*(pi*S_N2_n+(1-pi)*S_N2_d);

%-------------------------------------------------------------------------
% asset market clearing
%-------------------------------------------------------------------------
F(16) = - A + (1-e2)*eC*dA     + e2*dA2;
F(17) = - B + (1-e2)*(1-eC)*dB + e2*dB2;

%-------------------------------------------------------------------------
% liquidity premia
%-------------------------------------------------------------------------
F(18) = L_A - L_A2;
F(19) = L_B_rho_bar - L_B2_rho_bar;

%=========================================================================
% subfunctions
%-------------------------------------------------------------------------
function u = u(q,sigma)
% CRRA utility: u(q)
if sigma < 1
    u = q.^(1-sigma)./(1-sigma);
elseif sigma == 1
    u = log(q);
end

function mu = mu(q,sigma)
% marginal utility: u'(q)
mu = q.^(-sigma);

function w = w(q,sigma,theta)
% w function in effective bargaining power
w = theta + (1-theta)*mu(q,sigma);
